%% Stock Return Extrapolation, Option Prices, and Variance Risk Premium
% By Adem Atmaz
% This file generates Table 3: Predictive power of variance risk premium
% It also requires 1 function file: olshac.m
% Default code generates Univariate and Joint regression results
% Typically it takes less than 3 minutes to run

clear all;  clc;
format compact
format

% Parameter Values follow from  "Atmaz_Table_4_ParameterValues.m"
alpha=0.50;
theta=0.25;
mu=0.0212;
kappa=0.3567;
zeta=0.0044;
sigma=0.1625;
rho=0;
r=0.0056;

D0=100;
a=2.75;
vec_theta=[0, 0.25, 0.50, 0.75];
N_th=length(vec_theta);

% Simulations and Regressions
tau=1/12;  % variance horizon in years
h=3;       % predictive regression holding return horizon in months 1,3,12
% Simulation quantities
N_sim=1000;              % ##        % number of simulations: use 1000
S_end_year=1000;                     % sample size
int = 1/252;           %             % time interval for each year in a path
time_grid = 0:int:S_end_year;        % Time vector from 0 to end year
N_time=length(time_grid);
N_rets=N_time-1;                    % Number of returns in each path
H=1/12;     % Data frequency in the regression e.g. 1/12 for monthly
N_pp=H/int; % number of data points for each data point in regression: 21 for month

% omitting the first burn out period
keep_final_years=S_end_year; % keep all of it
Obs=keep_final_years*12;

% 1 path vectors
vec_D = zeros(1,N_time); % Allocate output vector
vec_V = zeros(1,N_time); % Allocate output vector
vec_X = zeros(1,N_time); % Allocate output vector

% All paths monthly matrices
N_month=(12*S_end_year)+1;
SM_V = zeros(1,N_month); % Allocate output vector
SM_X = zeros(1,N_month); % Allocate output vector
SM_D = zeros(1,N_month); % Allocate output vector
SM_S = zeros(1,N_month); % Allocate output vector
SM_VSR = zeros(1,N_month); % Allocate output vector
SM_VRP = zeros(1,N_month); % Allocate output vector
SM_RV  = zeros(1,N_month); % Allocate output vector

% Univariate regression of VRP
vec_beta_VRP        =zeros(N_sim,N_th);
vec_tstat_hac_VRP  =zeros(N_sim,N_th);
vec_R2_adj_VRP      =zeros(N_sim,N_th);

% Joint regression of VRP and RV
Joint_vec_beta_VRP    =zeros(N_sim,N_th-1);
Joint_vec_beta_RV     =zeros(N_sim,N_th-1);
Joint_vec_tstat_VRP_hac =zeros(N_sim,N_th);
Joint_vec_tstat_RV_hac  =zeros(N_sim,N_th);
Joint_vec_R2_adj        =zeros(N_sim,N_th-1);


% Begin simulation and regression
rng('default'); % fix the random number generator  seed   for replication
vec_Z1=randn(N_sim,N_time);
vec_Z2=randn(N_sim,N_time);
tic
for s = 1:N_sim
    s
    for th=1:N_th
        theta=vec_theta(th);
        % Stock price quantities
        M_DELTA=((kappa+(rho*sigma)*(1+theta))^2)-((2+theta)*(1+theta)*sigma^2);
        M_c=theta/(alpha*(1+theta));
        M_b=-(2+theta)/(kappa+((rho*sigma)*(1+theta))+sqrt(M_DELTA));
        M_a=a;
        M_sigma_1=(1+(rho*sigma*M_b))/(1-alpha*M_c);
        M_sigma_2=(sqrt(1-rho^2)*sigma*M_b)/(1-alpha*M_c);
        M_Var_con=(M_sigma_1^2)+(M_sigma_2^2);
        M_Cov_con=(rho+(sigma*M_b))/(1-alpha*M_c);
        % VSR and VRP quantities
        M_zeta_v=zeta*M_Var_con;
        M_kappa_v=kappa+(sigma*M_Cov_con);
        M_eta_v=(theta/M_b)*(1-(alpha*M_c))*M_Var_con;
        
        v0=M_zeta_v/kappa; % LR stock variance level
        M_RN_LRM_vt=(M_zeta_v+(r*M_eta_v))/(M_kappa_v+(0.5*M_eta_v));
        
        tt1=M_kappa_v+alpha;
        kappa1=(0.5*tt1)-(0.5*sqrt((tt1^2)-(4*alpha*(M_kappa_v+(0.5*M_eta_v)))));
        kappa2=(0.5*tt1)+(0.5*sqrt((tt1^2)-(4*alpha*(M_kappa_v+(0.5*M_eta_v)))));
        EE0=(1-exp(-kappa*tau))/(kappa*tau);
        EE1=(1-exp(-kappa1*tau))/(kappa1*tau);
        EE2=(1-exp(-kappa2*tau))/(kappa2*tau);
        
        % RV A and B
        RV_A=v0*(1-EE0);
        RV_B=EE0;
        
        % VSR A and B and C
        VSR_A=(M_RN_LRM_vt*(1-EE1-((kappa1/(kappa2-kappa1))*(EE1-EE2))))...
            +((M_zeta_v/(kappa2-kappa1))*(EE1-EE2));
        VSR_B=EE1-(((M_kappa_v-kappa1)/(kappa2-kappa1))*(EE1-EE2));
        VSR_C=(M_eta_v/(kappa2-kappa1))*(EE1-EE2);
        
        % VRP A and B and C
        VRP_A=(v0*(1-EE0))-VSR_A;
        VRP_B=EE0-VSR_B;
        VRP_C=-VSR_C;
        
        % Initalization of each path
        D0=100;
        V0=zeta/kappa; % long run mean of  fundamental  variance
        X0=(r+(0.5*M_Var_con*V0))/(1+theta); %  inital  X is equal to its long run mean        
        vec_D(1)=D0;
        vec_V(1)=V0;
        vec_X(1)=X0;
        
        for j = 1:N_rets
            Z1=vec_Z1(s,j);
            Z2=vec_Z2(s,j);
            vec_V(j+1) = max(0, vec_V(j))+((zeta-(kappa*max(0, vec_V(j))))*int)...
                +(sigma*rho*sqrt(max(0, vec_V(j)))*sqrt(int)*Z1)...
                +(sigma*sqrt(1-rho^2)*sqrt(max(0, vec_V(j)))*sqrt(int)*Z2);
            vec_X(j+1) = vec_X(j)+((r+(0.5*M_Var_con*vec_V(j))-(1+theta)*vec_X(j))*alpha*int)...
                +(alpha*M_sigma_1*sqrt(max(0, vec_V(j)))*sqrt(int)*Z1)...
                +(alpha*M_sigma_2*sqrt(max(0, vec_V(j)))*sqrt(int)*Z2);
            vec_D(j+1) = vec_D(j)+((mu*vec_D(j))*int)+(vec_D(j)*sqrt(max(0, vec_V(j)))*sqrt(int)*Z1);
            
        end
        
        % Get end of month values for each process from daily values
        SM_V= vec_V(1:N_pp:end);
        SM_X= vec_X(1:N_pp:end);
        SM_D= vec_D(1:N_pp:end);
        SM_S= SM_D.*exp(M_a+(M_b*SM_V)+(M_c*SM_X)); %stock price
        
        % Simulated VSR VRP and RV
        SM_RV=(RV_A+(RV_B*M_Var_con.*SM_V))*(10000/12); %normalized to get monthly values in percentage squared
        SM_VSR=(VSR_A+(VSR_B*M_Var_con.*SM_V)+(VSR_C.*SM_X))*(10000/12);  %normalized to get monthly values in percentage squared
        SM_VRP=SM_RV-SM_VSR;        
        
        % Keep only the selected final years and omit the first years as burnout if needed
        SM_RV=SM_RV(end-Obs:end);
        SM_VRP=SM_VRP(end-Obs:end);
        SM_S=SM_S(end-Obs:end);
        SM_D=SM_D(end-Obs:end);
        
        % Get h period ahead log returns
        SM_Ret_h=log(SM_S(1+h:end))-log(SM_S(1:end-h)); % log ret
        %SM_Ret_h=(SM_S(1+h:end)./SM_S(1:end-h))-1;      % returns
        SM_Ret_h_ann=(12/h)*SM_Ret_h*100;  %annualized in percentage terms        
        
        % Univariate Regression of h period returns on VRP
        % ===============================================        
        reg_Y=SM_Ret_h_ann';
        LT=length(reg_Y);
        reg_X=[ones(LT,1) SM_VRP(1:end-h)'];
        lag=round(0.75*(LT^(1/3)));    %  recommended hac lag choice
        
        [beta, R2, R2adj, X2_NW, X2_HH, X2_R, std_NW, std_HH, std_R, t_NW, t_HH, t_R]=olshac(reg_Y,reg_X,lag,lag);
        vec_beta_VRP(s,th)=beta(2);
        vec_tstat_hac_VRP(s,th)=t_NW(2);
        vec_R2_adj_VRP(s,th) = R2adj;        
        
        % Joint Regression of h period returns on VRP and RV
        % ===============================================
        if theta > 0   % joint regression except for benchmark model theta=0
            reg_X=[ones(LT,1) SM_VRP(1:LT)' SM_RV(1:LT)'];
            [beta, R2, R2adj, X2_NW, X2_HH, X2_R, std_NW, std_HH, std_R, t_NW, t_HH, t_R]=olshac(reg_Y,reg_X,lag,lag);
            Joint_vec_beta_VRP(s,th)=beta(2);
            Joint_vec_beta_RV(s,th)=beta(3);
            Joint_vec_tstat_VRP_hac(s,th)=t_NW(2);
            Joint_vec_tstat_RV_hac(s,th)=t_NW(3);
            Joint_vec_R2_adj(s,th) = R2adj;            
        end  % end if for joint regression
        
    end   % end theta
    
end    % end sim
clc;
toc

%  REPORT   UNIVARIATE REGRESSION on VRP QUANTITIES
tstat_hac_VRP=round(mean(vec_tstat_hac_VRP),2)
R2_VRP=round(mean(vec_R2_adj_VRP)*100,2)

%  REPORT   JOINT REGRESSION QUANTITIES
Joint_tstat_VRP_hac=round(mean(Joint_vec_tstat_VRP_hac),2)
Joint_tstat_RV_hac=round(mean(Joint_vec_tstat_RV_hac),2)
R2_Joint=round(mean(Joint_vec_R2_adj)*100,2)

clear alpha D0 EE0 EE1 EE2 EEv j jj kappa kappa1 kappa2 m M_a M_b M_c;
clear M_corr_sv M_Cov_con M_DELTA M_eta_v M_kappa_v M_Ob_LRM_vt M_r M_RN_LRM_vt;
clear M_sigma_1 M_sigma_2 M_sigma_v  M_Var_con M_zeta_v mu n rho sigma;
clear theta tt1 v V0 Z1 Z2 zeta;
clear vec_Z1 vec_Z2;

%% THE END










